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ABSTRACT 


The mechanical and rheological properties of an MR fluid depend on the induced 
microstructure of the imbedded ferrous particles. When subject to an externally applied 
magnetic field these particles magnetize and align themselves in chains parallel to the 
applied field. The microstructure of these chains is a function of several parameters 
including particle size, applied magnetic field strength, and viscosity and velocity of the 
surrounding fluid. This thesis will create a model from a first principle approach to 
accurately predict the microstructure in a variety of different situations. The model 
investigated assumes the particles become magnetic dipoles upon the application of the 
magnetic field and that particle interaction is due solely to dipole-dipole interaction. Due 
to the inherently small size of the particles, drag is modeled using Stokes’ drag. This 
mathematical model will be used to create a computer simulation to visualize and analyze 
the subsequent transient microstructures formed. The model will assume a constant 
magnetic field applied (i.e., no spatial or time gradients) and that the effects of this field 


are felt instantaneously. 
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I. INTRODUCTION 


A. CONTROLLABLE FLUIDS 


Magnetorheological (MR) and electrorheological (ER) fluids are a class of 
“smart” materials that are characterized by their ability to reversibly transform from 
liquid state to a Bingham solid. They are fluids that have either magnetically permeable 
(or electrically conductive) microscopic particles suspended in them. The transformation 
from liquid state to Bingham solid occurs by the application of a magnetic (or electric) 
field to the fluid. This magnetic (or electric) field causes the suspended particles to align 
in chains along the field lines in a manner to reduce the overall energy of the field. The 
existence of these chains changes many bulk properties of the fluid. Of practical interest 


is the change in viscosity which causes the fluid to behave like a Bingham solid. 
The Bingham model used for modeling MR fluids relates the total shear stress 7 
to the shear rate y and H (magnitude of applied magnetic field) according to the equation 


r=| real 








ery 


where 7,(H) is a yield stress that is a function of the applied magnetic field and 77 is the 


composite bulk viscosity of the fluid [1]. This equation is phenomenological in nature 
where the values in the equation are determined experimentally instead of being deduced 


from a first principle approach. 


Recently it has been found that a more detailed approach to predicting the 
behavior of MR fluids has become necessary due to the limitations of the above approach 
[1]. First, the Bingham model is a macro scale approach (the fluid and particles are 
treated as a single continuum instead of a composite system) with no differentiation with 
particle level. The coupling between mechanical behavior and the magnetic field takes 
place at the particle level and is governed by first principles (conservation of momentum, 


Maxwell’s Laws, etc.). The Bingham model then is limited to a narrow range of 


1 


applicability commensurate with the experimental data. Second, the Bingham model 
tends to be inaccurate at low value of stress. In current applications where MR devices 
are used as feedback controls, low value stresses are important and the above model 
proves unsatisfactory. Third, the Bingham model is only applicable to 1-D simple shear 
flows with a transverse magnetic field applied. This is inadequate for multi-degree of 
freedom MR devices that are currently being designed. These reasons encourage a 


different model to be developed that is based on first principles. 


B. INDUSTRIAL APPLICATIONS 


The American inventor, Willis Winslow, was the first to recognize how to create 
a smart fluid using these principles [2, 3]. He did much of the initial pioneering work on 
ER fluids in the 1940s. Later, Jacob Rabinow investigated the same phenomenon using a 
magnetic field for use in a magnetic field clutch [2] and is considered the first to develop 
MR fluids. Although their works were conducted over half a century ago, it has only 
been recently that the use of these smart fluids has become more common in industrial 
applications. This is due primarily to the stability and durability requirements of modern 


designs [3]. 


Today the uses for these smart, controllable fluids are numerous and varied. One 
primary use is in hydraulic dampers and brakes. Because of the ability to rapidly change 
the working fluid viscosity, one has the ability to change the damping coefficient in 
dampers to give much better dynamic response and control. Other applications include 
better feedback to control items such as joysticks, responsive personnel armor, and MR 


polishing machines [4]. 


C. TYPICAL MR FLUID ARRANGEMENT 


A typical arrangement for an industrial application of a MR fluid is shown below. 


Electromagnet 


MR Fluid 


Electromagnet 


Field Lines 





Figure 1. MR Fluid Device Arrangement 


The MR fluid, consisting of a carrier fluid (usually a silicone oil) and the 
suspended particles (typically fine ferrous particles), is sandwiched in a small gap 
between two electromagnets. These magnets, when energized, create a magnetic field 
perpendicular to the flow of the MR fluid which causes the imbedded particles to form 
chains parallel to the applied field. The dynamic response of these particles in both a 


static fluid and a moving fluid is investigated in this thesis. 
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I. MAGNETIC FORCE 


A. DIPOLE MODEL 


When the magnetic field is applied to the MR fluid, the ferromagnetic particles 
become magnetized and interact with the field and with each other. The exact solution to 
these interactions is difficult (if not impossible) and involves the integration of Maxwell’s 
stress tensors across the entire volume of the magnetized solution. Instead, some 
simplifying assumptions need to be made in order to allow an easily calculable analytical 


solution without sacrificing accuracy. 


The first assumption to be made is to use a dipole model for the magnetized 
particles. This assumes that the particles are dumbbell in shape, with length L. One end 
contains the positive (North) pole and the other contains the negative (South) pole of the 
magnet as shown in Figure 2. The magnitude of the pole strength is denoted by g and 


arises because the applied magnetic field magnetizes the particle. 














Figure 2. Dumbbell Shape 


When the above dumbbell is placed in a magnetic field H, it experiences a torque 
about its center as described by the formulaz = LqgH sin(@), where @is the angle between 
the magnetic field and dumbbell. The quantity L*g is given a special name, the magnetic 
moment, and is denoted by m. The model used in this thesis uses the dipole model and is 
defined by determining the value of m in the limit where L goes to zero but the torque 
remains finite. This is the case of magnetic spheres which are used in most MR 


applications. 


The magnitude of the dipole moment determines the interaction of the particles 
with the external magnetic field and the interaction of the particles with each other. It is a 
function of the magnitude of the applied external field (H with units Amp/m), the volume 
of the particles, and the magnetic permeability of both the particles and the surrounding 


fluid [5]. Specifically it is given by the relation 


m= (4ma°) as ad He (1) 
Hy + 2fy 


where a = radius of the particle (meters), “4, = magnetic permeability of the particle 
(henry/meter), 44 = magnetic permeability of the fluid, and H),. = magnetic field at dipole 


location (amp/meter). 


Several assumptions need to be made when using the above formula. The 
presence of a dipole alters the magnetic field in its vicinity (this is what causes particles 
to interact with each other) and this implies that the value of H),. needs to be calculated at 
every point. However, this variation is assumed to be negligible when calculating the 
magnetic dipole since the external fields applied are relatively large (on the order of 200 
kA/m) and the variations caused by the dipoles are several orders of magnitude smaller. 
Therefore Hj- is assumed to be equal to the applied magnetic field. The second 
assumption concerns the value of “4. Because the particles are ferrous, 4, is not a 
constant but varies with the applied magnetic field. However since the range of the 
applied field is small, often a fixed value, an average value of yu is used based on the 


values of the applied field. 
B. FORCE ON A DIPOLE 


The force on a dipole in a magnetic field is given by the product of the dipole 


moment and the gradient of the magnetic field as given by the expression 


4 OH 
Fo =m—— 
r Ai (2) 


a. . ree OH. eT eee: 
where F’ is the force in some arbitrary direction r and ae is the spatial derivative in the 
Tr 


r direction. One method that presents itself in determining the forces is to simply 
calculate the magnetic field at every point. Theoretically this could be done by 
calculating the external applied field and modifying it by the perturbations caused by the 
presence of the dipoles. Since the location of the dipoles constantly changes, this 
calculation would have to be performed at every time step. Using this calculation (which 
would have to be performed numerically) the gradient at every point could be calculated 
and then the force on every particle could be determined. In reality this calculation is 
difficult to perform, requires specific algorithms for determining the field and the 
gradients, and requires massive computing power. A more simplistic approach was used 


for this model. 


The first assumption for a more simplistic approach is that the applied magnetic 
field is uniform. This assumption is valid since the applied magnetic field is enacted 
rapidly (assumed instantaneous), does not vary with time, and the fluid gap is very small 
compared to the surface area over which the field is enacted. Since the magnetic field is 
assumed uniform in space and time there is no gradient and the particles experience no 
net force due to the external field. The only magnetic force the particles experience is 


due to their mutual interactions. 


Consider two magnetic dipoles of identical strength at arbitrary positions rand 


Bs A magnetic field H, is applied parallel to the z axis. The force between the two 





dipoles is given by 
7, = #t1 (1-3008" g, \z,—sin (2, )ay] 5 


if 


where f, is the force on particle 7 from particle j, 7, = i, =r, 


i 





, 6, is the angle from the z 


is a unit vector parallel to r,, and é, is a unit vector parallel to 


axis and 7, € iy? 


r 


e. “ (e, x A, [6]. This is shown in the below figure. 


oi 





Figure 3. Relationship Between Two Particles 


This equation models the interaction between dipoles and it is useful to examine 
this equation quantitatively to obtain a feel for the dynamics of the particles. By 
combining equations (1) and (2), it becomes apparent that the interaction force is 
proportional to the square of the applied field, proportional to the square of the particles 
volume, and the direction of the force is a function of the relative location of the two 
particles. This last item is what causes the particles to form stable chains when the 


magnetic field is applied. Examine only the radial term in equation (3). If@, is less than 
~54.6 degrees, the particles tend to attract. Otherwise they tend to repel. 


It is more useful to transform equation (3) into Cartesian coordinates. Using the 


same x, y, and z directions as shown in Figure (2) and defining Q = 3m’ u > the x, y, and z 


components are given as 





2 
it~ A 2 15 (4) 





jee 
a r |r. ) 
Y Y yy 
2 
Q SZi | Zi 
fg a| Oe (6) 
ij i i 


where fjjx is the x component of fj, fijy is the y component of fj, and fj, is the z component 


of f;;. The derivations for the above equations are attached in Appendix A. 


In a suspension of N particles each with an assumed induced magnetic 
dipole moment of m, the total magnetic force due to dipole interaction on a particle i is 
the sum of the contributions of all of the other particles in the suspension. In algebraic 


form 


N 
Fux, = > Fie (7) 


i#j 
N 
Fiuiy = DS Siy (8) 
iF] 
N 
Fy. = 2 Sis (9) 
iF] 


where F,,,,is the total magnetic force on particle i in the x direction, Fj, 1s the total 
magnetic force on particle i in the y direction, and F,,,,is the total magnetic force on 


particle i in the z direction. To determine the dynamic behavior of the particles in the 
fluid, these equation are calculated at every time step, the deviation in the current position 
of the particles are calculated, the values of the forces are recomputed at the next time 
step based on the particles’ new position, and the process is repeated until the end of the 


computational time. 
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Hi. OTHER FORCES 


A. DRAG FORCE 


The drag force on a spherical particle moving in a viscous fluid is a function of 
the pressure difference across the sphere (form drag) and the surface shear stress (viscous 
drag). In general this expression can be complicated to solve. In the specific case of the 
small particles used in MR fluids a number of simplifying assumptions can be made to 
more easily determine this drag force. The flow can be assumed to be laminar due to the 
small clearances between the electromagnets which the fluid flows between and the small 
velocities analyzed in this thesis. Another simplifying case arises due to the small 
particle size (~micro meter) which implies that the Reynolds number based on diameter 
(Rep) is less than 1. Both of these assumptions allow the viscous drag force to be 


modeled by the well understood Stokes’ drag which is given by 
F.=62na a (10) 
dt 


where F’ is the drag force in the r direction, 7 is the viscosity of the fluid and a is the 


radius of the particle. 
B. GRAVITY 


The force of gravity is neglected in this model based on the fact that the 
gravitational force that would tend to make the particles settle is a much weaker force 
than the magnetic force that acts between the dipoles. This is obviously not true when no 
magnetic field is applied, but the gravitational settling is ignored by assuming that the 
suspension is thoroughly mixed before the application of the field and that the particles 


are randomly distributed in the carrier fluid. 
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C. BROWNIAN MOTION 


Brownian motion is characterized by the random walk of particles in a fluid due 
to the bombardment of molecules. In determining if Brownian motion should be 
considered in any model of MR fluids, this effect should be compared to other effects 
which determine the dynamic behavior of the particles. When there is no applied 
magnetic field this is not the case, however by assuming that the fluid is mixed 
immediately before the field is applied would negate any effects of Brownian motion. 
Once the field is applied the relative effect of Brownian motion compared with the 
magnetic forces can be determined by analyzing the coupling constant / which is defined 
as the ratio of the interaction energy of two dipoles in contact and the thermal energy [6]. 
Specifically 


_ EF Moya’ HH” 
a a 9k,T 





(1) 


where yz, is the magnetic permeability in a vacuum, a is the particle radius, H is the 
magnetic field, yis the magnetic susceptibility of the particle, k,is the Boltzmann 


constant and T is the absolute temperature. In all cases modeled in this thesis 2 >>1 and 


Brownian motion is ignored. 
D. REPULSIVE FORCES 


The particles themselves and any walls that physically constrain the MR fluid are 
modeled as hard surfaces. Therefore a fictitious repulsive force must be modeled to 
ensure that a particle in physical contact with either a wall or another particle behaves as 
hard. The characteristics of this force are such that when two particles touch (the 
distance between two dipoles is 2a apart) the repulsive force exactly balances the 
attractive force between the dipoles, and when the distance between the dipoles is greater 


than 2a the force is negligibly small. The proposed form of this force is given below 


| ee Pane (12) 
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where K, and K2 are constants to be determined. The exponential term was chosen to 
give a function that rapidly decays as rj increases and is a commonly used mathematical 


model for these types of interactions [7]. 


To determine Kj, apply the condition that when rj is equal to 2a then the repulsive 
force must equal the attractive force between the dipoles. From equation (3) the dipoles 


attract when 6, is less than ~54.6 degrees and the attractive force also causes 6, to tend to 


zero. This is what causes the particles to align in chains that characterize the MR fluid. 


Assume that the particles will touch when @, is small. From Figure 2, this implies that 
x,, and y, are negligibly small. Appling this condition to equations (4-6) gives 


523 Ja, 
fix =9, fiy =O, and f= 9)2-“B |, 


ij I, 


When the particles are touching z, =7; =2a which, when combined with the 


above equation, gives 


Zo 
(2a)’ 





Fiz =i 


Combining this result with equation (9) gives a value for K, = —. 
a 


The constant K> is determined by a much less rigorous means. It must be negative 
to give a decaying characteristic and its magnitude is selected by a trial and error 
approach. On one hand a high magnitude gives a steeper decay which is advantageous 
since this more closely approximates reality. However, if the value is too large, the 
repulsive term can become extremely large for small distances and leads to numerical 


instabilities. A value of K, =—12 was chosen as a balance between these two competing 


factors based on numerical experiments. This gives a steep decay while allowing a more 


manageable time step. 
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Using the values of K; and K>2 and transforming equation (9) into Cartesian 
coordinates gives the following expressions for the repulsive force on a particle i from 


particle j in the x, y, and z directions as 





ee = Od, (13) 
8a i; 
QO 1 2] y; 

Teesty = o4 e : (14) 
8a , 
Q 1 ] Zi 

Tessctes = 8q* e ee (15) 


For the physical interaction with the walls of the container, a similar approach 
was taken, and the equation developed for the interaction of a particle with the 
floor/ceiling is given as 


30 2, -30| (4-4) 
Sz -= e “ti » | “ | (16) 





where /,;,is the force on particle i from the wall in the z direction, z, is the absolute 


distance from the bottom boundary to particle 7 in the z direction, and H is the total height 
of the volume (not to be confused with the use of H elsewhere as the magnetic field 
strength). Equations identical to (16) are used for the horizontal boundaries with the 


substitutions for the particles x and y positions (instead of z,) and the length and width of 


the containment area (instead of H). 
E. INTERACTION WITH ELECTROMAGNET 


Up to this point the discussion of the physics of the interactions of the particles in 


a MR fluid is exactly the same as if it was an ER fluid (replace the electromagnet with 
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charged parallel plate conductors and replace some of the magnetic constants with their 
electrical equivalents). A primary difference between the two arises in the physics of the 
interaction between the electromagnet (MR) and the interaction with the charged 
conducting plate (ER). In the latter case, the charged plates induce an electric dipole 
(exactly analogous to the electromagnets inducing a magnetic dipole), but the electric 
dipoles interact with the plates in an easily understood manner. The presence of the 
electric dipoles themselves will induce a current distribution on the plates and then these 
dipoles are electrically attracted to the plates because of this current distribution. A well 
documented manner to calculate this interaction is by the method of image charge [5]. 
Basically the interaction of a dipole that is a distance L from the plate is identical to 
assuming there is an infinite number of equal dipoles on the other side of the plate at 
distances nL where n=1,2,3..... Therefore an electric dipole will interact with the 
conducting plate, specifically will be attracted to the plate and attach itself to the plate. If 
there are multiple dipoles, they will form chains, and the chains will anchor themselves to 
the plate and behave as if the chain was infinitely long. This is what allows an ER to 


have a shear stress; the chains are anchored to the plate. 


There is no analogy in the MR case. There is no such thing as a magnetic current 
produced at the boundary of the electromagnet that would allow the use of the dipole 
image method to determine the interaction of the chain with the magnet [6]. Another way 
to look at this is to consider a single magnetic dipole between the magnets. Assuming a 
constant magnetic field, the dipole would not be attracted to either of the magnets. There 
seems to be nothing to lock the chain in place and therefore an MR fluid would not be 
able to have a shear stress. Experimentally, this is not true. The chains do become 


locked. 


There are two reasonable theories as to how this locking occurs. The first is to 
question the assumption that the field is uniform. Away from the edges of the magnets, 
due to the small gap between the magnets, it is safe to assume a uniform field. When the 
magnets are close together the field lines away from the edges do not spread out and 
consequently there is no gradient. However, at the edges of the magnets, this is not true. 
The fields bulge outward and tend to wrap around, causing large gradients. One proposal 
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is that away from the edges of the magnets, the chains are free to move (are not locked to the 
magnet), but as the bulk fluid flow sweeps them toward the edge of the magnet, they become 
locked in this area of high field gradients and effectively form a lattice type wall. Other 


chains being swept along will then build up behind this lattice wall. 


Another theory approach is to again to question the uniformity of the field, this time 
at the fluid/magnet interface. To explain this effect requires the use of Maxwell’s laws. The 


equation of specific use is 
V-B=0 (17) 
where B is magnetic flux density (Tesla). The relationship between H and B is given by 
B=wyH (18) 


where 4 is the magnetic permeability of the substance through which H exists. Using 


equations (17) to solve for the normal component of B across the discontinuity between the 


magnet and the fluid (there are no tangential components) implies 
n(,H, — 44H,)=0 or 44H, = 4A, 


where the subscripts refer to the magnetic permeability and magnetic field of the magnet and 
fluid respectively [8]. This shows at the interface between the magnet and fluid there is a 
jump discontinuity in the magnetic field (assuming that the magnetic permeability of the two 


materials are not equal). 


The model presented here assumes the second explanation for a physical interaction 
between the magnet and dipoles. The exact force caused by this gradient is unknown, but it 
is assumed that it is incredibly short ranged, and that causes a force of attraction such that, 
when multiplied by the frictional coefficient between the particle and magnet, leads to a 
frictional force that is substantially larger in magnitude as compared to the force that tends to 
sweep the particle along. This assumption is valid since the force tending to sweep the 
particle along with the flow is a function the flow velocity at the particle’s location. Since 
the particle is small (~5 microns) and resides at the interface, the flow velocity is 
approximately zero (no slip condition). In other word, a particle that happens to touch the 
magnet becomes locked in place, but particles in the stream, away from the wall, do not 
experience this force. 
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IV. MODEL FOR INTERACTION 


A. NEWTON’S SECOND LAW 


The description of a particle’s motion in a MR fluid can be determined using 
Newton’s second law of motion. In formulating a model for the motion of an arbitrary 


particle i apply this law in the x direction as follows 


d°x, 
DF =" (19) 





where the left hand side of the equation is the sum of all of the forces on particle i in the x 
direction and m is the particle’s mass (not dipole moment). The left hand side includes 
the dipole interaction force, drag force, and repulsive forces due to contact with other 
particles and the walls. Combining equations (7), (10), (13) and (16) gives the following 


second order differential equation to solve 


d°x. | 
m—7 +D— =F, (20) 





N 
where D = 67a, and F,, = Fy, + 2 He etek Dre 
i#j 
Equation (20) can be solved numerically, in its current form, using a range of 
techniques (for instance a Runge-Kutta algorithm). To make the computations more 
simple, integrate equation (20) over a sufficiently small time step 7 such that the term on 
the right hand side can be assumed constant. This gives an equation for the change in the 


position of the particle in the x direction during this time step 7 as shown below 


ee Voit -= 
Dae D l-e (21) 





Ly 


where V, is the velocity of the particle in the x direction at the beginning of the time step. 
This equations shows that if 7 is chosen such that it is several orders of magnitude less 


than ue then the second term on the right hand side can be ignored and therefore 











Ft 
Ax _ 1X 
D (22) 
and similarly 
A vy 
= D (23) 
jeder 
Az=— 
D (24) 


for the y and z components. For the MR fluids analyzed here a typical value of D is 
about 5E-7 sec’' which makes the time step on the order of 1E-9 sec. In reality this will 
be the upper limit on the time step. Initially a much smaller time step will be used in the 
computer simulation. This is due to how the program randomly establishes the initial 
positions of the particles and that they tend to overlap. A smaller time step is required to 


“push” the particles off of each other and the wall without destabilizing the computations 


with excessively large positional changes at each time step. 
B. STATIC FLUID MODEL 


A computational algorithm was written to determine the dynamic motion of the 
particles in a MR fluid. The program takes user inputs for the length, width and height of 
the MR fluid area, the number of particles to simulate, the magnitude of the applied 
magnetic field (program assumes the direction is in the negative z direction), and the 
number of time steps to perform the algorithm. The program then randomly distributes 
the particles inside of the fluid area. Using this distribution the initial spacing between all 


of the particles is computed (the values of x, , y,,,z, and7,). Then the value of the dipole 


strength and drag coefficient is computed. Using the spacing between particles and the 
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FE 


value of the dipole moments, F, aa 


AS and F.. are calculated for every particle. If the 
particle is near the upper or lower wall (a distance between a and 1.3a) the forces are 
assumed zero for the reason of the interaction with the magnet/fluid boundary discussed 


in Chapter II]. Then the forces are computed using equations (20). A time step is 


computed based on the value of ‘D 2s described in the previous section. The updated 


position of every particle is then calculated using equations (22-24) and this new position 
is stored and plotted graphically if desired. This process is repeated for every time step 


using the updated positions from the previous time step. 


As discussed above, a minor issue arises in the initial random spacing, especially 
at higher particle densities. Sometimes the particles are randomly placed such that two or 
more particles are spaced where the distance between them is less than their diameter 
length apart or such that the spacing between a particle and a wall is less than the 
particles radius length apart. This is not physically possible since the particles and the 
wall are hard. To overcome this, the time step chosen for the first 10 time steps is several 
orders of magnitude less than what is used for the remainder of the computation. This 
allows the repulsive force terms to “push” the particles away from each other and the 
wall without creating an abnormally large positional change that would eject them from 
the MR fluid domain. A copy of the computer code used is attached in Appendix B with 


a more detailed discussion as to the inner workings. 
C. DYNAMIC FLUID MODEL 


The programs constructed to compute the dynamic motion of particles in a 
dynamic fluid are very similar to the one for the static case with a few alterations. Two 
separate programs were created, one for pressure driven flow and the other for shear 


driven flow (these were the only two specific flow types analyzed). 


In the pressure driven case (parallel flow with a parabolic velocity distribution) it 
is assumed that the flow velocity is in the x direction, does not vary in the x and y 
directions and varies with a parabolic distribution in the z direction. The user inputs the 


meanline (maximum) flow and the program computes the value of the velocity at every 
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point in the MR fluid. Using this velocity distribution, another term is added to the right 
hand side of equation (20) to account for the drag force due to the fluid flow. Equation 
(20) now becomes 

d°x, dx, 


Me ae (25) 








where U;is the flow velocity at the position of particle 7. Note that since the flow is only 
in the x direction no modification is required for the equations of motion in the y and z 
directions. Applying the same arguments above that allowed for the inertial term to be 
ignored allows for the computation of the change in the position of the particle in the 


same manner as for the static case. 


The program for the shear driven flow (Couette flow) is the same as for the 
pressure driven flow, but here the user specifies the flow velocity at the upper plate. The 
program then calculates the velocity at all other points in the fluid and simulates the 


particle motion exactly the same as for the pressure driven flow. 
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V. SIMULATION RESULTS 


A. STATIC FLUID SIMULATION 


The first qualitative analysis to examine is the effect on time response of the fluid 
as a function of particle density. The results shown are for three simulations where all 
parameters are held constant with the exception of particle density. The various 
parameters used are shown in the below table. In all cases the size of the rectangular 


volume is 100 X 100 X 100 micrometers with hard walls bounding the area. The 


magnetic field is applied in the negative z direction. The value of D used to calculate 


the time step has a value of 1.744E-7 based on the below parameters. 

















Number of Fluid Viscosity Fluid Particle Applied 

Particles Permeability Permeability Magnetic Field 
Simulation 1 40 .25 Pas 1.26E-6 N/A’ .00377 N/A? 200 kA/m 
Simulation 2 70 .25 Pas 1.26E-6 N/A’ | .00377 N/A’ 200 kA/m 
Simulation 3 100 25 Pas 1.26E-6 N/A’ | .00377 N/A? 200 kA/m 




















Table 1. Parameters for Simulations 1, 2 and 3 


The simulations were conducted out for 100,000 time steps which equates to a 
simulation time of 1.7 milliseconds. The figures below show the particle microstructures 


at various times for the above simulations in both a 3-D and top down view. 
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Simulation 1, Initial Particle Distribution, 3-D View 














Figure 4. 
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Simulation 1, Initial Particle Distribution, Top Down View 
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Figure 5. 














Simulation 1, Time = 0.16 milliseconds, 3-D View 


Figure 6. 

















Simulation 1, Time = 0.16 milliseconds, Top Down View 


Figure 7. 
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Simulation 1, Time = 1.7 milliseconds, 3-D View 


Figure 8. 

















T T T 
I I I 
I \ \ 
\ I I 
I I \ 
\ I al o 
rn en 2 ee nee \ Ca) le es 5 eae 
I I | 
\ \ a 
\ I \ 
\ I item 
Piss cape tle Co cep de ec il Sa Ml SUNS Sc Ne, IS oe ITS s_ _\ | 0 
I \ {I ) |o 
} i he 
I I i 
\ I \ 
I I \ N 
 ecetiaetiaths Sitietiaeth 4-554 
\ I I ad 
\ I I 
\ I \ 
\ I Na 
I I I o 
pe ie ee rs nan io) 
I I I 
\ I I 
i id i 
\ | I 
Peo Sew ee hee eI eel eee SS as 2a) WO 
\ \ \ oO 
a \ \ 
I I \ 
te I I 
I are \ + 
nee eat Came fi ys ae $= = 6S 
\ I \ I I \ 
\ I I I \ 
\ \ I I \ 
I \ \ I \ 
ice ce clece sete cose li eeme tel coe Mees ol ono sg Toesge che. esse coe, | 88 
| \ \ (=) 
I I I 
\ I \ 
I I \ 

I I Je a 
eee ae ia ie ie ae ee ee ae (“y1s 
bo J 
I I lia 
\ gor \ 

I I I - 
Pe Ee ee Ne eS PO 1 V6 

| = \ 

\ I | 

I I I 

I I I 

L { L fo) 
= oe co Az ad 

oO io) io) 


x 10° 


Simulation 1, Time = 1.7 milliseconds, Top Down View 


Figure 9. 
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Simulation 2, Initial Particle Distribution, 3-D View 


Figure 10. 
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Simulation 2, Initial Particle Distribution, Top Down View 


Figure 11. 
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Simulation 2, Time = 0.16 milliseconds, 3-D View 


Figure 12. 

















Simulation 2, Time = 0.16 milliseconds, Top Down View 


Figure 13. 
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Simulation 2, Time = 0.86 milliseconds, 3-D View 


Figure 14. 
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Simulation 2, Time = 0.86 milliseconds, Top Down View 


Figure 15. 
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Simulation 2, Time = 1.7 milliseconds, 3-D View 


Figure 16. 
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Simulation 2, Time = 1.7 milliseconds, Top Down View 


Figure 17. 
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Simulation 3, Initial Particle Distribution, 3-D View 


Figure 18. 

















Simulation 3, Initial Particle Distribution, Top Down View 


Figure 19. 
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Simulation 3, Time = 0.16 milliseconds, 3-D View 


Figure 20. 
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Simulation 3, Time = 0.16 milliseconds, Top Down View 


Figure 21. 


30 

















Simulation 3, Time = 0.86 milliseconds, 3-D View 


Figure 22. 
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Simulation 3, Time = 0.86 milliseconds, Top Down View 


Figure 23. 
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Simulation 3, Time = 1.7 milliseconds, 3-D View 


Figure 24. 

















Simulation 3, Time = 1.7 milliseconds, Top Down View 


Figure 25. 
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Qualitatively, from the above figures, it is apparent that the speed at which the 
particles form chains is a function of the density of the particles in the fluid. Examining 
the two extreme cases (Simulations 1 and 3), much longer structures with fewer smaller 
structures are evident in Figure 20 than exist in Figure 6. Therefore, this model 
demonstrates the experimentally verified fact that particle density is a factor in response 


time of the MR fluid [3]. 


The second qualitative analysis to examine is the effect on time response of the 
fluid as a function of applied magnetic field strength. The results shown are for three 
simulations where all parameters are held constant with the exception of the magnetic 
field. The various parameters used are shown in Table 2. In all cases the size of the 


rectangular volume is 100 X 100 X 100 micrometers with hard walls bounding the area. 


The magnetic field is applied in the negative z direction. The value of D used to 


calculate the time step has a value of 1.744E-7 based on the below parameters. 



































Number of Fluid Fluid Particle maine 
Particles Viscosity Permeability Permeability Field 
Simulation 1 70 .25 Pas 1.26E-6 N/A? | .00377 N/A? 150 kA/m 
Simulation 2 70 .25 Pas 1.26E-6 N/A* | .00377 N/A? 200 kA/m 
Simulation 3 70 .25Pas 1.26E-6 N/A* | .00377 N/A? 250 kA/m 
Table 2. Parameters for Simulations 4, 5 and 6 
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Simulation 4, Initial Particle Distribution, 3-D View 


Figure 26. 

















Simulation 4, Initial Particle Distribution, Top Down View 


Figure 27. 
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Simulation 4, Time = 0.86 milliseconds, 3-D View 


Figure 28. 

















Simulation 4, Time = 0.86 milliseconds, Top Down View 


Figure 29. 
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Simulation 4, Time = 1.7 milliseconds, 3-D View 
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Simulation 4, Time = 1.7 milliseconds, Top Down View 
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Figure 31. 























Simulation 5, Initial Particle Distribution, 3-D View 


Figure 32. 
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Simulation 5, Initial Particle Distribution, Top Down View 


Figure 33. 
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Simulation 5, Time = 0.86 milliseconds, 3-D View 


Figure 34. 
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Simulation 5, Time = 0.86 milliseconds, Top Down View 


Figure 35. 
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Simulation 5, Time = 1.7 milliseconds, 3-D View 


Figure 36. 

















Simulation 5, Time = 1.7 milliseconds, 3-D View 


Figure 37. 
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Simulation 6, Initial Particle Distribution, 3-D View 





Figure 38. 
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Simulation 6, Initial Particle Distribution, Top Down View 


Figure 39. 
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Simulation 6, Time = 0.86 milliseconds, 3-D View 


Figure 40. 
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Simulation 6, Time = 0.86 milliseconds, Top Down View 


Figure 41. 
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Simulation 6, Time = 1.7 milliseconds, 3-D view 


Figure 42. 
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Simulation 6, Time = 1.7 milliseconds, Top Down View 











Figure 43. 
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Qualitatively, from the above figures, it is apparent that the speed in which the 
particles form chains is a function of the strength of the applied magnetic field. 
Examining the two extreme cases (Simulations 4 and 6) it is clear that the structures are 
much closer to completion in Figure 40 than those in Figure 28 (both are 0.86 
milliseconds into the simulation). This result is intuitive based on the fact that the 
strength of the dipole force is proportional to the square of the applied magnetic field 


strength (Equations | and 3). 
B. DYNAMIC FLUID SIMULATION 


The dynamic fluid simulation is where the real use for the MR model is realized. 
Most MR applications in industry use the MR effect on a moving fluid. It is desirable for 
the dynamic model to be able to accurately predict the microstructures of the particles in 
the MR fluid and from there predict the apparent viscosity and shear stress of the MR 
fluid. One theory for predicting the shear stress has already been developed and makes 
its predictions based on the chain density in the MR fluid and the angles the chains make 


in relation to the moving fluid [9]. 


A simulation of a MR fluid in shear was run in order to show that these 
measurements could be made in order to test the model against experimental evidence 
and to use the model to help design MR fluid devices. The simulation results are shown 


below and the following parameters were used. 











at Fluid Fluid Particle Me Velocity of 
Pare Viscosity Permeability | Permeability Field Top Plate 
Simulation 7 40 .25Pas | 1.26E-6 N/A® | .00377 N/A* | 250 kA/m 0.1 m/s 























Table 3. Parameters for Simulation 7 


The following figures show the distribution of particles and their orientation in a 
dynamic flow. The dimensions for the area were changed in order to create a higher 


particle density without adding more particles to the volume. In this case the size of the 
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rectangular volume is 50 X 100 X 50 micrometers with the longer direction in the 
direction of the fluid velocity. The particles were only seeded in the left half of the 


volume to allow the fluid to push the particles to the right without interference from the 


wall. The magnetic field is applied in the negative z direction. The value of ” D used to 


calculate the time step has a value of 1.744E-7 as was the case for the above simulations. 
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Figure 44. Side View of a MR Fluid in Shear 


In the above figure the fluid is flowing from the left to the right based on the shear 
force developed due to the top plate moving to the left at 0.1 m/s. The angle that the 
chains make with relation to the fluid vary, but could easily be measured and averaged. 
Examining the top down view of the fluid shown below, would allow for the 


determination of the chain density. 
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Figure 45. Top View of the MR Fluid in Shear 


The dynamic model allows the user to easily input common variable parameters 
in a program in order to determine the dynamic microstructure in two common flow 
conditions (linear flow from shear and parabolic flow from a pressure gradient). For 
other, more permanent parameters (particle or fluid magnetic permeability, fluid 
viscosity, etc.) individual lines of code which set these parameters must be altered. For a 
more detailed description of what parameters the user inputs when running the code and 


other parameters directly set in the code see Appendix B. 


The pressure driven, parabolic velocity profile fluid simulation is shown below. 
Unfortunately it is harder to determine the microstructure, in this case than for the shear 
case. The red line inserted on the below figures shows one chain and how it bulges in the 
middle based on the higher flow velocity there. The parameters for this simulation are 
exactly the same as those shown in Table 3, except the flow velocity is the centerline 


flow, not the flow at the top wall. 
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Side View of MR Fluid with Parabolic Flow 


Figure 46. 
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Top Down View of MR Fluid with Parabolic Flow 


Figure 47. 
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VI. CONCLUSION 


The goal of the model developed was for it to be simple, easy to use, require little 
empirical data (based on first principles) and be accurate. The model satisfies the 
requirement to be simple. It uses very well understood laws (Newton’s Second Law, 
dipole interaction force, and Stokes’ drag) to describe the forces and accelerations of the 
particles. Through a mathematical justification it ignores the inertial mass of the particles 
(being dominated by the viscous forces) to simplify the calculations even further. The 
physical interaction between the particles and the walls was chosen to be in a form that 
would balance out the other dominant forces in a way that was short ranged. This 


technique is also commonly used in other similar types of models. 


The programs were written in order allow multiple simulations with a minimum 
amount of work. Instead of the user having to laboriously input every parameter required 
for every simulation the programs only require the user to input a small number of 
parameters that were assumed to be the most varied. The disadvantage to this approach is 
that the user must modify the code in order to change parameters such as particle or fluid 
magnetic permeability, fluid viscosity, or particle size. However, it was assumed that 
these parameters would not be changed as frequently as the magnetic field strength, size 
of the volume, or number of particles, so they were not requested by the program for 


every simulation. 


The model used is based entirely on first principles so no empirical data is 
required for the simulation of particles. This could change if the model requires 
modification in order to accurately predict experimental results. For instance, in some 
regimes, the Stokes’ approximation may no longer be valid. In this case the model 
should be modified to more accurately describe the shear and pressure forces on a 
submerged body and this may have to be done with empirical data. It is the hope and 


belief that this will not be required, but it is a possibility. 


As far as accuracy is concerned the best way to check the model is against 


carefully controlled laboratory experiments. Qualitatively the model matches observed 


47 


data such as the chain formation and the time scale in which these structures form 10. 
However, there are several assumptions made that may have to be modified. The first 
assumption that limits the applicability of the dynamic model is the characterization of the 
velocity profile. In reality, the formation of particle chains alters the imposed flow. This is 
how MR fluids are able to withstand a shear and alter the bulk viscosity of the fluid. The 
model presented does not allow for the changing of the velocity profile. This is not a simple 
problem to solve, because the only way to determine how the flow changes based on the 
particle formations and how the particle formations vary based on the flow is to numerically 
solve to sets of equations simultaneously. Either Euler’s equations or the Navier-Stokes’ 
equations must be solved at each time step along with the other equations to determine the 
forces acting on the particles. This would require the integration of the model for the 
dynamic behavior of the magnetic particles with a numerical solver for fluid dynamic. The 
location of the particles at each time step would represent a boundary in the flow that would 
be solved with a flow solver. The model presented here does not attempt to perform any type 
of alteration of the flow based on the particle dynamics. Therefore, the model is not good for 
long simulation times, but it is still valid in the short term (before the initial velocity of the 
fluid is altered). In other words, this model accurately describes the initial particle dynamics, 
but is poor in the limit where the initial flow velocity would have been modified by the 


particle structures. 


Another area requiring more detailed study is in the interaction between the particle 
chains and the magnet providing the magnetic field. As discussed earlier, this interaction is 
well understood in the ER case, but not for the MR. Models for ER fluid accurately describe, 
based on the interaction with the chains and the electrode, how particle chains will merge to 
form even larger structures around one second after the field is applied. These larger 
structures are also observed in MR fluids, but the method in which they form is being 
debated. Since it is doubted that the chains interact with the magnet in the same way that the 
ER fluid’s chains interact with the electrodes, the same process is not occurring. Some have 
proposed that Brownian motion needs to be included. While the Brownian motion has much 
smaller interaction energy than magnetism, it is postulated that Brownian interaction could 
cause the chains to bulge and this temporary, minor change in position of the chain could 
allow for nearby chains to attract this chain. Again, this is a postulation, and more study is 


required for the particle-magnet interaction. 


48 


APPENDIX A. DERIVATIONS OF EQUATION 4, 5, AND 6 


Based on the figure below define ra as the component of iE in the radial 


direction and f; as the component of i. in the angular direction as shown. 


y 





Figure 48. Geometrical Relationship Between Two Particles 


After dropping the subscripts for convenience and decomposing /f. into its 


r 


Cartesian components gives 
f, = f, sin(@) sin(¢) 


f. = f, sin() cos(g) 
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f. = f, cos(@) 


where f,, f. and f, are the magnitudes of the components of the force in the y, x and z 
directions respectively. Similar decompositions of f, into its Cartesian components 


gives 
Ff, = fy c0s(@) sin(g) 
Ff, = fg c0s(O) cos(¢) 
f, =—f, sin(@) 
with the same definitions as above. Adding the components together gives 
f, = f, sin(@) sin(g) + f, cos(@) sin(g) 
f. = f, sin(@) cos(g) + f, cos(A) cos() 
f, = f, cos(@) — f, sin(@) . 


Substituting the definition of f, and f, from equation (3) and the geometric identities 


2 2: 
RE) ee 


7 
sin(@) => — 
x+y 
cos(@) = i 
. 
cos(~) = 


Xx 


into the above equation obtains 


J -2)125 2 


r r r 
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APPENDIX B. COMPUTER CODE AND DESCRIPTION 


Below is the computer code for the computation of the static flow problem and a 


description of the various sections. 


sThesis program 

clear all 

a=5*10%-6; %m radius of particle 
Vol=pi*a%*3*4/3; SVolume of particle 








mass=7850*Vol; %Skg mass of particle 

tf=input ("Number of time steps ub 

height=input ('Height of volume (micro meter)'); 
height=height*10%*%-6; %converts to meters 

length=input ('Length of volume (micro meter)'); 
length=length*10%-6; %Sconverts to meters 

width=input ('Width of volume (micro meter)'); 
width=width*10%-6; 

N=input ('Number of particles'); 

H=input ('Magnetic Field intensity (kA/m)'); %~200 kA/m 





H=H*1000; 


xXinit=length*rand(N,1); Sinitial x dist of particles 


yinit=width*rand(N,1); 
Zinit=height*rand(N,1); 


Sinitial 


y dist of particles 





Sinitial z dist of particles 








vis=.25; Sfluid viscosity [Pa*s 

uf = 1.257E-6; %Spermeability of fluid 

up = .00377; %permeability of particle 

m = (4/3) *pi*H*a%3* (up-uf) / (up+2*uf); Smagnetic moment 


D=6*pi*vis*a; 
tcheck=mass/D; 


Stokes drag force coefficient 
Sintrinsic time scale 
































Q = 3*m*2*uf; 

ysys = zeros(tf,N); ty position of particles, column 1 refers to 
particle 1, column 2 refers to particle 2, ect 

xsys = zeros(tf,N); %x position of particles, column 1 refers to 
particle 1, column 2 refers to particle 2, ect 

zsys = zeros(tf,N); %z position of particles, column 1 refers to 
particle 1, column 2 refers to particle 2, ect 

delx = zeros(N,N); %Sdifference in the x position between the particles 
dely = zeros(N,N); %Sdifference in the y position between the particles 
delz = zeros(N,N); sdifference in the z position between the particles 





r=zeros (N,N); 
Fmx=zeros (N,N 


Fpx=zeros 
Fpy=zeros 
Fpz=zeros 
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Fmz=zeros (N,N); 





“5* (5* (delz(i,j)/xr(i,j)) 





(1, 3))) *exp (- 


i from particle j in x 


ysys(1,:) = yinit'; 
xsys(1,:) = xinit'; 
zsys(1,:) = zinit'; 
for t = 1:tf 
xX = xsys(t,:); 
y = ysys(t,:); 
Z = zsys(t, ; 
for i=1:N 
for jJ=l:N 
delx(i,j) = x(i)-x(9); 
dely(i,j) = y(i)-y(3); 
delz(i,j) = z(i)-z(9); 
(i,j) = sqrt (delx(i,j)*2+dely (i,j) *2+delz(i,4j)%*2) 
if i==j 
Fmx (1,3) =0; 
Fmy (i, 3) =0; 
Fpy (i, 3)=0; 
Fpx (i, j)=0; 
Fpz(i,4)=0; 
Fmz (i,j) =0; 
else 
Fmx (i, 3) =-Q*delx(i,j)/r ( 
1); Sforce on particle i from eee i 
force 
Fpx (i, j)=2*Q/ ((2*a) *4) * (delx (i,j) /(r 
12* ((r (i, j)/(2*a))-1)); %Sforce on particle 
direction due to physical interaction meray 






































Fmy (i, 3) =-O*dely (i,j) /r ( “5* (5* (delz (i,j) /r(i,j))*2- 
1); Sforce on particle i from particle j in i direction due to magnetic 
forces 
Fpy (i, 3) =2*Q/ ( (2*a) *4) * (dely (i, 3) / (r(i, 4) )) *exp (- 
12* ((r(i,j)/(2*a))-1)); Sforce on particle i from particle j in y 
direction due to physical interaction (collision) 
Fmz (i, 3)=-Q*delz (i,j) /r(i, 4) *5* (5* (delz (i, 3) /zr (i,j) )*%2- 
3); 
Fpz (i,j) =2*O/ ((2*a) *4) * (delz (i,j) /(r (4,34) )) *exp (- 
12* ((x(i,j)/(2*a))-1)); 
Fx=Fmxt+tFpx; %Stotal force in x direction 
Fy=FmytFpy; %stotal force in y direction 
Fz=FmztFpz; %Stotal force in z direction 
end 
end 
Ftx=sum(Fx,2); Stotal force on each particle in x direction 
Fty=sum(Fy,2); stotal force on each particle in y direction 
Ftz=sum(Fz,2); %Stotal force on each particle in z direction 
end 
for gq=1:N 
Ftz (q) =Ftz (q)+2*Q/ ( (2*a) *4) *exp (-30* (z(q) /(2*a)-.5) ); 
Ftz (q) =Ftz (q) -2*Q/ ((2*a) *4) *exp (-30* ( (height-—z (q)) / (2*a) 
-5)); Sloop incorperates force due to repulsion of wall 
Ftx (q) =Ftx (q) +2*Q/ ( (2*a) *4) *exp (-30* (x(q) / (2*a)-.5)); 
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A2- 
in x direction due to magnetic 


Ftx (q) =Ftx (q) -2*Q/ ( (2*a) *4) *exp (-30* ( (length-x (q) ) / (2*a) - 


5)); 
Fty (q) =Fty (q) +2*Q/ ( (2*a) *4) *exp (-30* (y(q) / (2*a)-.5)); 
Fty (q) =Fty (q) -2*Q/ ((2*a) *4) *exp (-30* ( (width-y (q)) / (2*a) - 
5)); 
end 
Lf t< LQ 


tau=tcheck/1E7; 
elseif t<100 
tau=tcheck/1E4; 
elseif t<1000 
tau=tcheck/1E2; 




















else 
tau=tcheck/10; 
end 
xsys(tt1,:) Ftx'*tau/D; 














x 
ysys(tt1,:)=ytFty'*tau/D; 
zsys(t+1,:)=z+Ftz'*tau/D; 








plot3(xsys(t+1,:),ysys(tt+1,:),zsys(tt1,:),'0', 'Markersize',25); 


axis ([{0,length,0,width,0,height]) 
pause(.0001) 
end 


The first section clears all currently stored variables, sets particle size and 


calculates the mass of the particles. 


SThesis program 

clear all 

a=5*10%-6; %m radius of particle 
Vol=pi*a*3*4/3; %Volume of particle 
mass=7850*Vol; %Skg mass of particle 

















The next section is for the user to input the number of time steps, volume 
dimensions, number of particles, and magnetic field strength. It also converts these 


inputs into the proper units. 








tf=input ('Number of time steps Ls 

height=input ('Height of volume (micro meter)'); 
height=height*10%-6; %converts to meters 
length=input ('Length of volume (micro meter)'); 
length=length*10*%-6; %Sconverts to meters 
width=input ('Width of volume (micro meter)'); 


width=width*10%*-6; 
N=input ('Number of particles'); 
H=input ("Magnetic Field intensity (kA/m)'); %~200 kA/m 


as, 





H=H*1000; 


The next section of code randomly decides the initial positions of the 


particles and stores them in the variables shown. 


Xinit=length*rand(N,1); Sinitial x dist of particles 
yinit=width*rand(N,1); %Sinitial y dist of particles 
Zinit=height*rand(N,1); Sinitial z dist of particles 
































The next section set some parameters of the MR fluid (magnetic 
permeability of the particle and fluid and fluid viscosity) and computes the magnetic 
dipole moment, Stokes’ drag coefficient, and the intrinsic time scale (to be used in 


determining the actual time scale later in the program). 


vis=.25; sfluid viscosity [Pa*s] 








uf = 1.257E-6; %permeability of fluid 
up = .00377; %Spermeability of particle 
m = (4/3) *pi*H*a%3* (up-uf) /(upt2*uf); tmagnetic moment 


D=6*pi*vis*a; SStokes drag force coefficient 
tcheck=mass/D; %intrinsic time scale 
Q = 3*m*2*uf; 


The next section preallocates memory for the matrices that will be used to 
either store or compute the motion of the particles. The matrices xsys, ysys, and zsys will 
store the actual positions of the particles. The columns refer to the individual particles 
(the first column is the position of particle 1, the second column refers to the position of 
particle 2, etc.) and the rows refer to the time (the first row is the initial distribution, the 
second row is the position of the particles after one time step, etc.). The matrices delx, 
dely and delx temporarily store the differences in the x, y and z direction between 
particles. The index of the matrix determines the difference in position of which 
particles. For example, delx(4,9) is the difference in the x position between particles 4 
and 9. These matrices get written over after each time step. The matrix r is similar to 
delx, dely and delz but stores the difference in the radial direction between particles. The 
matrices Fmx, Fmy, and Fmz store the x, y and z components of the force between two 
particles due to their magnetic dipoles. For example, Fmy(2,6) is the dipole force in the y 


direction between particles 2 and 6. Fpx, Fpy, and Fpz are similar to Fmx, Fmy, and Fimz 
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except that the former relate to the physical repulsive force due to the particles being hard 





spheres. 

ysys = zeros(tf,N); ty position of particles, column 1 refers to 
particle 1, column 2 refers to particle 2, ect 

xsys = zeros(tf,N); %x position of particles, column 1 refers to 
particle 1, column 2 refers to particle 2, ect 

zsys = zeros(tf,N); %z position of particles, column 1 refers to 
particle 1, column 2 refers to particle 2, ect 

delx = zeros(N,N); sdifference in the x position between the parti 
dely = zeros(N,N); *sdifference in the y position between the parti 
delz = zeros(N,N); %sdifference in the z position between the parti 





r=zeros (N, 


Fmx=Zeros 
Fmy=zeros 
Fmz=zeros 
Fpx=zeros 
Fpy=zeros 
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es 
es 
es 


Fpz=zeros 


The next section takes the initial particle distribution and stores this data in 


the first row of the corresponding position matrices. 


ysys(1,:) = yinit'; 
xsys(1,:) = xinit'; 
zsys(1,:) = zinit'; 


The next section is the heart of the program that calculates the dipole and 
physical interaction between particles. These nested loops calculate first the delx, dely, 
delz and r matrices for the current time step. Then set the diagonal elements of the force 
matrices to zero (a particle does not interact with itself). It then uses the equations 


discussed in the thesis to calculate all of the other elements in the force matrices. 





for t = 1:tf 
xX = xsys(t,:); 
y = ysys(t,:); 
Z= zsys(t,:); 
for i=1:N 
for j=1l:N 
delx(i,j) = x(1)-x(j); 
dely(i,j) = y(i)-y(3); 
delz(i,j) = z(1)-z(j); 
r(i,j) = sqrt (delx(i,4)*2+dely (i, 4)*2+delz(i,4)*2); 
if i==j 
Fmx (1, 3) =0; 
Fmy (i, 3) =0; 
Fpy (i, 3) =0; 


a 


Fpx (i, 4) =0; 
Fpz (i, 3) =0; 
Fmz (i,j) =0; 
else 
Fmx (i, J) =-Q*delx(i,j)/r ( “5* (5* (delz(i,j)/r(i,j))%*2- 


1); sforce on particle i from eee j in x direction due to magnetic 
force 














Fpx (i, j)=2*Q/ ( (2*a) *4) * (delx (i, 4) / (r (i, 3) )) *exp (- 
12* ((r(i,j)/(2*a))-1)); Sforce on particle i from particle j in x 
direction due to physical interaction (co ision) 

Fmy (i, 3) =-Q*dely(i,j)/r (i,j) *5* (5* (delz (i,j) /r (i,j) )%*2- 
1); Sforce on particle i from particle j in y direction due to magnetic 


force 











Fpy (i, 3) =2*Q/ ( (2*a) *4) * (dely (i, 3) / (x (i, 4) )) *exp (- 
12* ((r(i,j)/(2*a))-1)); Sforce on particle i from particle j in y 
direction due to physical interaction (collision) 

Fmz (i, 3) =-Q*delz (i, 4) /r(i,4)*5* (5* (delz (i, 3) /r (i, j))*2- 
3); 

Fpz (i,j) =2*O/ ((2*a) *4) * (delz (i,j) /(r (i, 4) )) *exp (- 


12* ((x (i, 3)/(2*a))-1)); 


























Fx=Fmxt+Fpx; %Stotal force in x direction 

Fy=FmytFpy; %total force in y direction 

Fz=FmztFpz; %Stotal force in z direction 

end 

end 
Ftx=sum(Fx,2); Stotal force on each particle in x direction 
Fty=sum(Fy,2); sStotal force on each particle in y direction 
Ftz=sum(Fz,2); stotal force on each particle in z direction 














end 


The next section incorporates a loop to add the force due to the physical 


interaction with the wall. 











for q=1:N 

Ftz (q) =Ftz (q) +2*O/ ( (2*a) *4) *exp (-30* (z (q) / (2*a i. ee 

Ftz (q) =Ftz (q) -2*Q/ ((2*a) *4) *exp (-30* ( (height~—z (q) ) / ( 45) ) 7 
Sloop incorperates force due a ss wie of wall 

Ftx (q) =Ftx (q)+2*Q/ ( (2*a) *4) *exp (-30* (x(q) /(2*a)-.5) ); 

Ftx (q) =Ftx (q) -2*Q/ ( (2* ie ) *exp (-30* ( (length-x (q) )/ (2*a)-.5)); 

Fty (q) =Fty (q) +2*Q/ ( (2*a) *4) *exp (-30* (y(q) / (2*a)-.5) ); 

Fty (q) =Fty (q) -2*Q/ ( (2*a) *4) *exp (—30* ( (width-y (q) )/(2*a)-.5)); 


end 


The next section uses the intrinsic time scale to create an actual time scale 
based on the time that the program is simulating. Extremely small time scales are used 
initially to allow the initial structures to begin forming and then the time scales are 
enlarged as the structures become closer to equilibrium. 
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SF -C<10 

tau=tcheck/1E7; 
elseif t<100 
tau=tcheck/1E4; 
elseif t<1000 
tau=tcheck/1E2; 














tau=tcheck/10; 
end 


The final section calculates the new position of every particle, stores it in 


the appropriate matrices and plots the positions on a 3-D graph. 








xsys(tt1,:)=x+Ftx'*tau/D; 
ysys(tt1,:)=ytFty'*tau/D; 
zsys(tt+1,:)=z+Ftz'*tau/D; 




















plot3 (xsys(t+1,:),ysys(tt+1,:),zsys(tt1,:),'0', 'Markersize',25); 


axis ([{0,length,0,width, 0,height]) 
pause (.0001) 


Below is the attached code for the shear flow problem. This code is very 
similar to the static code above, except for a few lines described. Sentences in red are the 


descriptions of the changes and do not appear in the code. 


sThesis program 

lear all 

a=5*10%-6; %m radius of particle 
Vol=pi*a%*3*4/3; SVolume of particle 
mass=7850*Vol; Skg mass of particle 





Q 





tf=input ('Number of time steps "); 

height=input ('Height of volume (micro meter)'); 
height=height*10%-6; %Sconverts to meters 
length=input ('Length of volume (micro meter)'); 
length=length*10%-6; %Sconverts to meters 
width=input ('Width of volume (micro meter)'); 
width=width*10%-6; 





U=input (‘Velocity of top plate /s)'); eee: 


N=input ('Number of particles'); 
H=input ("Magnetic Field intensity (kA/m)'); %~200 kA/m 
H=H*1000; 





Xinit=length*rand(N,1); tinitial x dist of particles 
yinit=width*rand(N,1); Sinitial y dist of particles 
zZinit=height*rand(N,1); 


ag 


vis=.25; 
uf = 1.257E-6; Spermeability of fluid 
up = .00377; spermeability of particle 








m = 4*pi*H*a%*3* (up-uf) / (upt+2*uf); Smagnetic moment, may to modify H 
since local magnetic field may not be applied field (see paper 
micromechanical model for MR fluids) 

D=6*pi*vis*a; SStokes drag force coefficient 

tcheck=mass/D; %intrinsic time scale 

Q = 3*m*2*uf; 
































ysys = zeros(tf,N); ty position of particles, column 1 refers to 
particle 1, column 2 refers to particle 2, ect 

xsys = zeros(tf,N); %x position of particles, column 1 refers to 
particle 1, column 2 refers to particle 2, ect 

zsys = zeros(tf,N); 

delx = zeros(N,N); 

dely = zeros(N,N); 

delz = zeros(N,N); 

r=zeros(N,N); 

Fmx=zeros (N,N); 

Fmy=zeros (N,N); 

Fpx=zeros (N,N); 

Fpy=zeros (N,N); 

Fpz=zeros (N,N); 

Fmz=zeros (N,N); 

ysys(1,:) = yinit'; 

xsys(1,:) = xinit'; 

zsys(1,?2) = zinit'; 


for t = 1:tf 
xX = xsys(t 
y = ysys(t,:); 
Z = zsys(t H 





delz (i 


if i==j 


a Se Lae tes es 


else 





Fmx (i, J) =-Q*delx(i,j)/r (i, 4) *5* (5* (delz (i, 4)/r (i,j) )%*2- 
1); Sforce on particle i from particle j in x direction due to magnetic 
force 
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Fpx (i, j) =Q/ ((2*a) *4) * (delx (i, 5) / ( (4,5) )) *exp (- 
12* ((r(i,j)/(2*a))-1)); Sforce on particle i from particle j in x 
direction due to physical interaction (collision) 


Fmy (i, 3) =-OQ*dely (i,j) /r(i, 4) *5* (5* (delz (i, 3) /xr(i,3)) *2- 


1); Sforce on particle i from particle j in y direction due to magnetic 
force 








Fpy (1,3) =Q/ ((2*a) *4) * (dely (i, 3) /(r (4,3) )) *exp (- 

12* ((r(i,j)/(2*a))-1)); Sforce on particle i from particle j in y 
direction due to physical interaction (collision) 

Fmz (i, 3) =-Q*delz(i,j)/r(i, 4) *5* (5* (delz (i,j) /r (i,j) )%*2- 








Sy 








Fpz (i,j) =2*O/ ((2*a) *4) * (delz (i,j) /(r (4,3) )) *exp (- 
12* ((r(i,j)/(2*a))-1)); 
Fx=Fmxt+Fpx; %Stotal force in x direction 
Fy=FmytFpy; %total force in y direction 
F2=Fmz+Fpz; 
end 








end 
Ftx=sum(Fx,2); Stotal force on each particle in x direction 
Fty=sum(Fy,2); stotal force on each particle in y direction 
Ftz=sum(Fz,2); %stotal force on each particle in z direction 
end 
for q=1:N 
if z(q)<1.2%*a && z(q)>=a 
Ftx(q)=0; 
Fty (q)=0; 
Ftz (q) =0; 
else 
Ftz (q) =Ftz (q) +2*Q/ ((2*a) *4) *exp (-30* (z(q) / (2*a)-.5 
Ftz (q) =Ftz (q) -2*Q/ ( (2*a) *4) *exp (-30* ( (height~—z (q) ) 
-5)); Sloop incorperates force due to repulsion of wall 
Ftx (q) =Ftx (q) +2*Q/ ((2*a) *4) *exp (-30* (x(q) / (2*a)-.5 
Ftx (q) =Ftx (q) -2*Q/ ((2*a) *4) *exp (-30* ( (length-x (q) ) 














Ftx (q) =Ftx (q) +D*U* (z (q) -a) /height; 









Fty (q) =Fty (q) +2*Q/ ((2*a) *4) *exp (-30* (y(q) / (2*a) -.5)); 
Fty (q) =Fty (q) -2*Q/ ((2*a) *4) *exp (-30* ( (width-y (q)) / (2*a) - 


end 
end 


if ¢<1Q 
tau=tcheck/1E7; 
elseif t<100 
tau=tcheck/1E5; 
elseif t<1000 
tau=tcheck/1E3; 








else 








tau=tcheck/100; 


end 
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xsys(tt1l,:)=x+Ftx'*tau/D; 
ysys(tt1,:)=ytFty'*tau/D; 
zsys(tt1l,:)=z+Ftz'*tau/D; 




















plot3 (xsys(t+1,:),ysys(tt+1,:),zsys(tt1,:),'0o', 'Markersize',25); 
axis ([{0,length,0,width,0,height]) 


pause (.0001) 
end 
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